####################################################################################
# Keeping Vigil: The Emergence of Vigilance Committees in Pre-Civil War America.  #
# Jonathan Obert & Eleonora Mattiacci                                             #
# Last modified: January 31, 2018                                                 #
###################################################################################

## This R file contains the code necessary to replicate the analysis in the main text.
## The analysis was carried out using R version 3.4.1 on a 3.5 GHz Intel Core i7 running OS Sierra 10.12.6


rm(list=ls())

#loading packages (function from Huff, Connor, and Joshua D. Kertzer. "How the Public Defines Terrorism." American Journal of Political Science 62.1 (2018): 55-71.)
ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
                        if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
                        sapply(pkg, require, character.only = TRUE)
}

packages <- c("foreign","survival", "simPH", "extrafont", "RColorBrewer", "MASS") 
ipak(packages)

#####
# Setting working directory
#####

setwd("") 
dat<-read.dta("data_replication.dta")
dat1<-subset(dat,dat$multiple==0);
dat1$slat<-dat1$sla_i*log(dat1$stop + 1)


#####
# Code to reproduce Model 1, Table 1
#####

model1m<-coxph(Surv(start, stop, committee)~  ef_i+ bre+ ef_i_bre +cluster(fips),data=dat1,method="efron");summary(model1m);cox.zph(model1m);
AIC(model1m);BIC(model1m);logLik(model1m)

#####
# Code to estimate effects:
#####

x<-(exp(model1m$coef[1])-1)*100;x
y<-(exp(model1m$coef[1]+model1m$coef[3])-1)*100;y

A<-((y-x)/y)*100;A

#####
# Code to reproduce Figure 1 in the Text
#####

lightgray<-"#d9d9d9";darkgray<-"#525252";gray<-"#666666"
col_seq<-c(as.vector(c(gray,darkgray,darkgray)),as.vector(c("black",lightgray,lightgray)))
line_seq<-c(as.vector(c(2,1,1)),as.vector(c(2,1,1)))

dat2m<-data.frame(ef_i=c(quantile(dat1$ef_i,probs=c(.25),na.rm=TRUE),quantile(dat1$ef_i,probs=c(.75),na.rm=TRUE)),bre=c(0,1),ef_i_bre=c((quantile(dat1$ef_i,probs=c(.25),na.rm=TRUE)*0),(quantile(dat1$ef_i,probs=c(.75),na.rm=TRUE)*1)))

plot(survfit(model1m,newdata=dat2m),ylim=c(0.85,1) ,conf.int=T,lty=c(as.vector(c(1,2,2)),as.vector(c(1,2,2))),lwd=line_seq,col= col_seq,ylab="Proportion of Counties With No Committee",xlab="Years till Committee Formation",xlim=c(0,11))
legend("bottomleft",c("No Social Frontier Conditions","95%C.I.","" ,"Social Frontier Conditions","95%CI"),lty=c(1,2,0,1,2),col=c(gray,darkgray,"","black",lightgray),lwd=c(3,3),border=lightgray,cex = .85,bty="o",box.col=lightgray)

######
# Code to reproduce Model 2, Table 1
#######

# We first ran the simplest model. A Ph test indicates the need to add an interaction of Slavery with Time
#model2<-coxph(Surv(start, stop, committee)~  ef_i+ bre+ef_i_bre +gini_i+manu_i+urb_i+sla_i+cluster(fips),data=dat1,method="efron");summary(model2);cox.zph(model2);

model2m<-coxph(Surv(start, stop, committee)~  ef_i+ bre+ef_i_bre +gini_i+manu_i+urb_i+sla_i+slat+cluster(fips),data=dat1,method="efron");summary(model2m);
AIC(model2m);BIC(model2m);logLik(model2m)

















